########Appendix Replication for Trussler BJPS 2019#########
rm(list=ls())
library(sandwich)
library(foreign)
library(stargazer)
library(car)
library(plm)
library(lfe)
library (sp)
library (rgdal)
library (raster)
library(maptools)
library(RColorBrewer)

#Change to working directory in which you have the data files
setwd("C:/Google Drive/Information and Accountability/Incumbency Bias Nationalization/Trussler BJPS Replication/ReplicationData")

##### Appendix 1 American National Election Study Data merged with District level Aggregate Data##### 

load("NES data for BJPS Appendix 1.Rdata")

#Models using district number of providers predicting individual level knowledge and attitudes

m.rate.ideo<-  felm(ideo.rate  ~ factor(year) + log(providers) +
                      factor(party.3) + age  + female + factor(race) +
                      factor(education) +  factor(income) + log(med.income) + log(pop) 
                    ,
                    data=nes, weights=nes$weight)
summary(m.rate.ideo)


m.rate.ideo.sy<-  felm(ideo.rate  ~  log(providers) +
                         factor(party.3) + age  + female + factor(race) +
                         factor(education) +  factor(income) +log(med.income) + log(pop)  | state.year |0 | state.year ,
                       data=nes, weights=nes$weight)
summary(m.rate.ideo.sy)


m.rate.ft<-  felm(ft.rate  ~ factor(year) + log(providers)  +
                    factor(party.3) + age + female + factor(race) +
                    factor(education) +  factor(income) + log(med.income) + log(pop)  ,
                  data=nes, weights=nes$weight)
summary(m.rate.ft)



m.rate.ft.sy<-  felm(ft.rate  ~  log(providers)  +
                       factor(party.3) + age + female + factor(race) +
                       factor(education) +factor(income) +log(med.income) + log(pop) | state.year | 0 | state.year,
                     data=nes, weights=nes$weight)
summary(m.rate.ft.sy)




m.know.house<-  felm(know.house  ~ factor(year) + log(providers)  +
                       factor(party.3) + age + female + factor(race) +
                       factor(education)+factor(income) +  log(med.income) + log(pop)  ,
                     data=nes, weights=nes$weight)
summary(m.know.house)

m.know.house.sy<-  felm(know.house  ~  log(providers)  +
                          factor(party.3) + age + female + factor(race) +
                          factor(education)+factor(income) +  log(med.income) + log(pop)  |state.year | 0 | state.year ,
                        data=nes, weights=nes$weight)
summary(m.know.house.sy)


m.know.senate<-  felm(know.senate  ~factor(year) +  log(providers)  +
                        factor(party.3) + age + female + factor(race) +
                        factor(education)+factor(income) +  log(med.income) + log(pop) ,
                      data=nes, weights=nes$weight)
summary(m.know.senate)

m.know.senate.sy<-  felm(know.senate  ~factor(year) +  log(providers) +
                           factor(party.3) + age + female + factor(race) +
                           factor(education)+factor(income) +  log(med.income) + log(pop) | state.year | 0 | state.year ,
                         data=nes, weights=nes$weight)
summary(m.know.senate.sy)

library(stargazer)

#Creating Table 1
stargazer(m.rate.ideo, m.rate.ideo.sy, m.rate.ft, m.rate.ft.sy,
          m.know.house, m.know.house.sy, m.know.senate, m.know.senate.sy,
          omit="year", star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          covariate.labels = c("log(Providers)","Independent", "Republican", 
                               "Age", "Woman", "Race-Other", "Race-White",
                               "Education - Grade School or Less", "Education - High School", 
                               "Education - Some College", "Income - 17-33 Percentile",
                               "Income - 23-67 Percentile", "Income - 68-95 Percentile",
                               "Income - 96-100 Percentile",
                               "log(Median Income)", "log(Population)"),
          dep.var.labels   = c("P(Rate Incumbent Ideology", "P(Rate Incumbent Thermometer)",
                               "P(Know House Majority)","P(Know Senate Majority)" ),
          model.numbers = F,
          add.lines = list(c("Year F.E.", "Yes","No", "Yes", "No", "Yes", "No", "Yes", "No"),
                           c("State-Year F.E.","No", "Yes", "No", "Yes", "No", "Yes", "No", "Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Political Knowledge",
          label="knowledge")


nes$split.ticket[nes$party.3=="Independent"] <- NA

m.split <-  felm(split.ticket  ~ factor(year) + log(providers)  + 
                   factor(party.3) + age + female + factor(race) + 
                   factor(education)+factor(income) +  log(med.income) + log(pop)  ,
                 data=nes, weights=nes$weight)
summary(m.split)



m.split.sy <-  felm(split.ticket  ~ log(providers)  + 
                      factor(party.3) + age + female + factor(race) +
                      factor(education)+factor(income) +  log(med.income) + log(pop) 
                    | state.year | 0 | state.year,
                    data=nes, weights=nes$weight)
summary(m.split.sy)



m.affect <-  felm(affective.polarization ~ factor(year) + log(providers)  + 
                    factor(party.3) + age + female + factor(race) +
                    factor(education)+factor(income) +   log(med.income) + log(pop) ,
                  data=nes, weights=nes$weight)
summary(m.affect)


m.affect.sy <-  felm(affective.polarization  ~ factor(year) + log(providers)    + 
                       factor(party.3) + age + female + factor(race) +
                       factor(education)+factor(income)+   log(med.income) + log(pop)  | state.year | 0 | state.year ,
                     data=nes, weights=nes$weight)
summary(m.affect.sy)


#Creating Table 2
stargazer(m.split,m.split.sy, m.affect, m.affect.sy, omit="year", star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          covariate.labels = c("log(Providers)","Independent", "Republican", 
                               "Age", "Woman", "Race-Other", "Race-White",
                               "Education - Grade School or Less", "Education - High School", 
                               "Education - Some College", "Income - 17-33 Percentile",
                               "Income - 23-67 Percentile", "Income - 68-95 Percentile",
                               "Income - 96-100 Percentile",
                               "log(Median Income)", "log(Population)"),
          dep.var.labels   = c("P(Split Ticket)", "Affective Polarization",
                               "P(Know House Majority)","P(Know Senate Majority)" ),
          model.numbers = F,
          add.lines = list(c("Year F.E.", "Yes","No", "Yes", "No" ),
                           c("State-Year F.E.","No", "Yes", "No", "Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Voting and Attitudes",
          label="polarization")


#Code to create Figure 1
coefficients <- c(coef(m.rate.ideo)["log(providers)"], coef(m.rate.ideo.sy)["log(providers)"], coef(m.rate.ft)["log(providers)"], coef(m.rate.ft.sy)["log(providers)"],
                  coef(m.know.house)["log(providers)"], coef(m.know.house.sy)["log(providers)"], coef(m.know.senate)["log(providers)"], coef(m.know.senate.sy)["log(providers)"],
                  coef(m.split)["log(providers)"],coef(m.split.sy)["log(providers)"], coef(m.affect)["log(providers)"], coef(m.affect.sy)["log(providers)"])

se <-c(sqrt(vcov(m.rate.ideo)["log(providers)","log(providers)"]), sqrt(vcov(m.rate.ideo.sy)["log(providers)","log(providers)"]), sqrt(vcov(m.rate.ft)["log(providers)","log(providers)"]), sqrt(vcov(m.rate.ft.sy)["log(providers)","log(providers)"]),
       sqrt(vcov(m.know.house)["log(providers)","log(providers)"]), sqrt(vcov(m.know.house.sy)["log(providers)","log(providers)"]), sqrt(vcov(m.know.senate)["log(providers)","log(providers)"]), sqrt(vcov(m.know.senate.sy)["log(providers)","log(providers)"]),
       sqrt(vcov(m.split)["log(providers)","log(providers)"]),sqrt(vcov(m.split.sy)["log(providers)","log(providers)"]), sqrt(vcov(m.affect)["log(providers)","log(providers)"]), sqrt(vcov(m.affect.sy)["log(providers)","log(providers)"]))

ci.90u <- coefficients + se*1.645
ci.90d <- coefficients - se*1.645
ci.95u <- coefficients + se*1.96
ci.95d <- coefficients - se*1.96


length(coefficients)
positions <- seq(-1,-12,-1)

positions <- c(-1,-2,-4,-5,-7,-8,-10,-11, -13,-14, -16,-17)


pdf(file="./AppendixFigures/IndividualLevelCoefficients.pdf", height=8, width=9)
par(mar=c(5.1,14,4.1,2.1))
plot(coefficients[seq(1,11,2)], positions[seq(1,11,2)], pch=16, xlim=c(-.25,.25), ylim=c(-18,1),
     axes=F, cex=1, xlab="Effect of 100% Change in Number of Broadband Providers", ylab="")
abline(v=0, lty=2, col="gray80")
abline(h=seq(-3,-15,-3), lty=2)

segments(coefficients, positions, ci.90u, positions , lwd=2)
segments(coefficients, positions, ci.90d, positions , lwd=2)
segments(coefficients, positions, ci.95u, positions , lwd=1)
segments(coefficients, positions, ci.95d, positions , lwd=1)
points(coefficients[seq(2,12,2)], positions[seq(2,12,2)], pch=17, cex=1)
axis(side=1, at=seq(-.2,.2,.1))
axis(side=2, at=seq(-1.5, -16.5 ,-3), las=2, label=c("P(Rate Incumbent Ideology)", "P(Rate Incumbent Thermometer)",
                                                     "P(Know House Majority)","P(Know Senate Majority)",
                                                     "P(Vote Split Ticket)","Affective Polarization"), col="white")
legend("topright", c("Year Fixed Effects", "State-Year Fixed Effects"), pch=c(16,17))    
dev.off()


##### Appendix 2 Broadband and Turnout####
rm(list=ls())
load(file="AggregateNatInc.Rdata")
load("Congressional Turnout Data.Rdata")

cd <- merge(cd, turnout, all.x=T, by.x="unique", by.y="district.years")


m.turnout <- plm(perc.turnout ~ factor(year)+log(providers) +log(med.income) + log(pop) ,
                 data=cd, model="within", index=c("cdfip","year"))
m.turnout.vcov <- vcovHC(m.turnout, type="HC0", cluster="group", method = "arellano")
m.turnout.se <- sqrt(diag(m.turnout.vcov))


#Creation of Table 3

stargazer(m.turnout, omit=c("cdfip", "year"),
          se=list(m.turnout.se),star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          covariate.labels = c("log(Providers$_{kt}$)", "log(Median Income$_{kt}$)", "log(Population$_{kt}$)"),
          dep.var.labels   = "% VAP Turnout$_{kt}$",
          model.numbers = F,
          add.lines = list(c("District F.E.", "Yes"),
                           c("Election Year F.E.","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Congressional Election Turnout",
          label="turnout")



m.incumbency1 <- plm(dv ~factor(year) + ptynow + dpres + inc3*log(providers) + 
                       log(med.income) + log(pop) +perc.turnout,
                     data=cd, model="within", index=c("cdfip","year"))
m.incumbency1.vcov <- vcovHC(m.incumbency1, type="HC0", cluster="group", method = "arellano")
m.incumbency1.se <- sqrt(diag(m.incumbency1.vcov))


m.incumbency2 <- plm(dv ~factor(year) + ptynow + dpres + inc3*log(providers) + 
                       log(med.income) + log(pop) +perc.turnout*log(providers),
                     data=cd, model="within", index=c("cdfip","year"))
m.incumbency2.vcov <- vcovHC(m.incumbency2, type="HC0", cluster="group", method = "arellano")
m.incumbency2.se <- sqrt(diag(m.incumbency2.vcov))


#Creation of Table 4

stargazer(m.incumbency1,m.incumbency2, omit=c("cdfip", "year"),
          se=list(m.incumbency1.se, m.incumbency2.se),star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          covariate.labels = c("Party Now$_{kt}$","% Democratic Vote for President$_{kt}$","Incumbency$_{kt}$",
                               "log(Broadband Providers$_{kt}$)","log(Median Income$_{kt}$)",
                               "log(Population$_{kt}$)", "% VAP Turnout$_{kt}$",
                               "Incumbency$_{kt}$*log(Broadband Providers$_{kt}$)",
                               "% VAP Turnout$_{kt}$*log(Broadband Providers$_{kt}$)"),
          dep.var.labels   = "% Democratic Vote for House$_{kt}$",
          model.numbers = F,
          add.lines = list(c("District F.E.", "Yes","Yes"),
                           c("Election Year F.E.","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Incumbency Advantage",
          label="incumbency")




m.incwin.pres1 <- plm(inc.vote ~ factor(year)*factor(party) +  log(providers)*pres.vote +
                        log(pop)  + log(med.income) +perc.turnout , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.pres1.vcov <-  vcovHC(m.incwin.pres1, type="HC0", cluster="group", method = "arellano")
m.incwin.pres1.se <- sqrt(diag(m.incwin.pres1.vcov))
summary(m.incwin.pres1)

m.incwin.pres2 <- plm(inc.vote ~ factor(year)*factor(party) +  log(providers)*pres.vote +
                        log(pop)  + log(med.income) +log(providers)*perc.turnout , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.pres2.vcov <-  vcovHC(m.incwin.pres2, type="HC0", cluster="group", method = "arellano")
m.incwin.pres2.se <- sqrt(diag(m.incwin.pres2.vcov))
summary(m.incwin.pres2)



m.incwin.party1 <- plm(inc.vote ~factor(year)*factor(party) + log(providers)*party.unity  + 
                         log(pop)  + log(med.income)+perc.turnout , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.party1.vcov <-  vcovHC(m.incwin.party1, type="HC0", cluster="group", method = "arellano")
m.incwin.party1.se <- sqrt(diag(m.incwin.party1.vcov))



m.incwin.party2 <- plm(inc.vote ~factor(year)*factor(party) + log(providers)*party.unity  + 
                         log(pop)  + log(med.income)+perc.turnout*log(providers) , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.party2.vcov <-  vcovHC(m.incwin.party2, type="HC0", cluster="group", method = "arellano")
m.incwin.party2.se <- sqrt(diag(m.incwin.party2.vcov))

#Creation of Table 5

stargazer(m.incwin.pres1, m.incwin.pres2,m.incwin.party1, m.incwin.party2, omit=c("icpsr.id","name","year","ptynow"),
          se=list(m.incwin.pres1.se, m.incwin.pres2.se,m.incwin.party1.se, m.incwin.party2.se),star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          order=c("pres.vote", "party.unity","log(providers)", "log(providers):pres.vote",
                  "log(providers):party.unity"),
          covariate.labels = c("Same Party Presidential Vote$_{jt}$",
                               "Same Party Presidential Vote$_{jt}$*log(Broadband Providers$_{jt}$)",
                               "Party Unity$_{jt}$",
                               "Party Unity$_{jt}$*log(Broadband Providers$_{jt}$)",
                               "log(Broadband Providers$_{jt}$)", 
                               "log(Population$_{jt}$)",
                               "log(Median Income$_{jt}$)",
                               "% VAP Turnout$_{jt}$",
                               "% VAP Turnout$_{jt}$*log(Broadband Providers$_{jt}$)"),
          dep.var.labels   = "% Incumbent Vote$_{jt}$",
          model.numbers = F,
          add.lines = list(c("Member F.E.", "Yes","Yes","Yes","Yes"),
                           c("Election Year*Party F.E.","Yes","Yes","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Incumbency Advantage",
          label="incvote")

##### Appendix 3 Strategic Politicians####
rm(list=ls())
load(file="AggregateNatInc.Rdata")

cd$uncontested <- NA
cd$uncontested[is.na(cd$inc.vote)]<- 1
cd$uncontested[!is.na(cd$inc.vote)]<- 0


keep <- c("unique","providers","cdfip","icpsr.id",
          "med.income","pop", "year")
cd.resid <- cd[keep]
cd.resid <-cd.resid[complete.cases(cd.resid),]

m.bb.inc <- plm(log(providers) ~ 0 + factor(year) + factor(icpsr.id) + 
                  log(med.income) +  log(pop) 
                , data=cd.resid, model="within", index=c("icpsr.id","year"))
m.bb.inc.se <- sqrt(diag(vcovHC(m.bb.inc,type="HC0",cluster="group",method="arellano")))

cd.resid$residual.inc<- residuals(m.bb.inc)

keep <- c("unique","residual.inc")
cd.resid <- cd.resid[keep]

cd <- merge(cd,cd.resid, all.x=T, by="unique")

m.resid.inc <- plm(residual.inc ~  factor(year) + factor(icpsr.id) +  inc.vote , data=cd,
                   model="within", index=c("icpsr.id","year"))
m.resid.inc.se <- sqrt(diag(vcovHC(m.resid.inc,type="HC0",cluster="group",method="arellano")))

m.resid.compete <- plm(residual.inc ~ factor(year) + uncontested, data=cd,
                       model="within",index=c("icpsr.id","year"))
m.resid.compete.se <- sqrt(diag(vcovHC(m.resid.compete,type="HC0",cluster="group",method="arellano")))

m.resid.pu <- plm(residual.inc ~ factor(year) + party.unity, data=cd,
                  model="within",index=c("icpsr.id","year"))
m.resid.pu.se <- sqrt(diag(vcovHC(m.resid.pu,type="HC0",cluster="group",method="arellano")))


stargazer(m.bb.inc, m.resid.inc,m.resid.compete,m.resid.pu,report = "vc*p")


#Creation of Table 6

stargazer(m.bb.inc, m.resid.inc,m.resid.compete,m.resid.pu, omit=c("cdfip","icpsr.id","year"),
          se=list(m.bb.inc.se, m.resid.inc.se,m.resid.compete.se,m.resid.pu.se),star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          covariate.labels = c("$ln(Med.Income_{it}$", "$ln(Pop_{it})$",
                               "$Inc.Vote_{it}$", "$Seat.Uncontested_{it}$","$Party.Unity_{it}$"),
          dep.var.labels = c("$ln(Providers_{it})$",
                             "$ln(Providers_{it}) -hat{ln(Providers_{it})}$"),
          model.numbers = F,
          add.lines = list(c("Incumbent F.E.", "Yes","Yes","Yes","Yes"),
                           c("Election Year F.E.","Yes","Yes","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Relation of Residual Variation to Electoral Marginality",
          label="residual")




##### Appendix 4 CPS Internet Supplement Validation #######
rm(list = ls())
load(file="CPS Internet Merged.Rdata")

cps.internet$providers_nonzero <- cps.internet$providers
cps.internet$providers_nonzero[cps.internet$providers_nonzero==0] <- 0.1
cps.internet$lnproviders <- log(cps.internet$providers_nonzero)


library(stargazer)
m.bb <- felm(broadband ~  factor(year)  +log(providers) +  lnpop+
               log(med_income)| county | 0 | county,
             data=cps.internet ,
             weights=cps.internet$wtsupp)

stargazer(m.bb, type="text")

#Creation of Table 7
stargazer(m.bb, omit=c("year"),
          star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          covariate.labels = c("$ln(Providers_{ict})$", "$ln(Population_{ict}$",
                               "$ln(Med.Income_{ict})$"),
          dep.var.labels   = "P(Home Broadband)",
          model.numbers = F,
          add.lines = list(c("County F.E.", "Yes"),
                           c("Year F.E.","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Broadband on Home Broadband Subscription",
          label="cps.broadband")



##### Appendix 5.1 Urban and Rural Districts ######

#Urban Rural districts seperately
rm(list=ls())
load(file="AggregateNatInc.Rdata")
urb <- read.csv("CensusUrbanRural.csv")

head(urb, 30)
names(urb)

urb$Geo_CD110[urb$Geo_CD110==0] <- 1
urb$cd <- urb$Geo_CD110
urb$state <- urb$Geo_STATE

urb$cd <- as.character(urb$cd)
urb$cd[nchar(urb$cd)==1] <- paste("0",urb$cd[nchar(urb$cd)==1], sep="" )
urb$state <- as.character(urb$state)
urb$state[nchar(urb$state)==1] <- paste("0",urb$state[nchar(urb$state)==1], sep="" )
urb$cdfip <- paste(urb$state,urb$cd, sep="")

keep <- c("cdfip","SE_T002_001","SE_T002_005")
urb <- urb[keep]

names(urb) <- c("cdfip","pop","rural")
urb$perc.rural <- urb$rural/urb$pop



keep <- c("cdfip","perc.rural")
urb <- urb[keep]

cd <- merge(cd, urb,by="cdfip", all.x=T)

#Define rural districts

rural.districts <- cd$cdfip[cd$year==2002 & cd$perc.rural>.50]


cd$rural[cd$cdfip %in% rural.districts]<-1
cd$rural[!(cd$cdfip %in% rural.districts)]<-0



m.incumbency2.urban <- plm(dv ~factor(year) + ptynow + dpres + inc3*log(providers)*rural + 
                             log(med.income) + log(pop) ,
                           data=cd, model="within", index=c("cdfip","year"))
m.incumbency2.urban.vcov <- vcovHC(m.incumbency2.urban, type="HC0", cluster="group", method = "arellano")
m.incumbency2.urban.se <- sqrt(diag(m.incumbency2.urban.vcov))
summary(m.incumbency2.urban)


#Creation of Table 8
stargazer(m.incumbency2.urban, omit=c("cdfip", "year"),
          se=list( m.incumbency2.urban.se),star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          covariate.labels = c("$Party_{jt}$","$D.Pres_{jt}$","$Incumbency_{jt}$",
                               "$ln(Providers_{jt})$","$ln(Med.Income_{jt})$",
                               "$ln(Population_{jt}$", 
                               "$Incumbency_{jt}*ln(Providers_{jt})$",
                               "$Incumbency_{jt}*Rural_j$",
                               "$ln(Providers_{jt})*Rural_j$",
                               "$Incumbency_{jt}*ln(Providers_{jt})*Rural_j$"),
          dep.var.labels   = "$D.House_{jt}$",
          model.numbers = F,
          add.lines = list(c("District F.E.", "Yes"),
                           c("Election Year F.E.","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Incumbency Advantage",
          label="urban.incumbency")

m.incwin.pres2.urban <- plm(inc.vote ~ factor(year)*factor(party) +  log(providers)*pres.vote*rural +
                              log(pop)  + log(med.income)  , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.pres2.urban.vcov <-  vcovHC(m.incwin.pres2.urban, type="HC0", cluster="group", method = "arellano")
m.incwin.pres2.urban.se <- sqrt(diag(m.incwin.pres2.urban.vcov))
summary(m.incwin.pres2.urban)

m.incwin.party2.urban <- plm(inc.vote ~factor(year)*factor(party) + log(providers)*party.unity*rural  + 
                               log(pop)  + log(med.income) , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.party2.urban.vcov <-  vcovHC(m.incwin.party2.urban, type="HC0", cluster="group", method = "arellano")
m.incwin.party2.urban.se <- sqrt(diag(m.incwin.party2.urban.vcov))



#Creation of Table 9
stargazer(m.incwin.pres2.urban, m.incwin.party2.urban, omit=c("cdfip", "year"),
          se=list(m.incwin.pres2.urban.se, m.incwin.party2.urban.se), star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          order=c("pres.vote", "party.unity","log(providers)", "log(providers):pres.vote",
                  "log(providers):party.unity"),
          covariate.labels = c("Same Party Presidential Vote$_{jt}$",
                               "Same Party Presidential Vote$_{jt}$*log(Broadband Providers$_{jt}$)",
                               "Same Party Presidential Vote$_{jt}$*Rural$_{j}$",
                               "Same Party Presidential Vote$_{jt}$*log(Broadband Providers$_{jt}$)*Rural$_{j}$",
                               
                               "Party Unity$_{jt}$",
                               "Party Unity$_{jt}$*log(Broadband Providers$_{jt}$)",
                               "Party Unity$_{jt}$*Rural$_{j}$",
                               "Party Unity$_{jt}$*log(Broadband Providers$_{jt}$)*Rural$_{j}$",
                               "log(Broadband Providers$_{jt}$)", 
                               "Rural$_{j}$)", 
                               "log(Population$_{jt}$)",
                               "log(Median Income$_{jt}$)",
                               "log(Broadband Providers$_{jt}$)*Rural$_j$" ),
          dep.var.labels   = "% Incumbent Vote$_{jt}$",
          model.numbers = F,
          add.lines = list(c("Member F.E.", "Yes","Yes"),
                           c("Election Year*Party F.E.","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Incumbency Advantage",
          label="urban.incvote")

##### Appendix 5.2 Cross-Sectional Models ######
rm(list=ls())
load(file="AggregateNatInc.Rdata")


#Incumbency Advantage, Cross Section Year


m.incumbency.02 <- felm(dv ~  inc3 + ptynow + dpres + inc3*log(providers) + 
                          log(med.income) + log(pop)|state| 0 | state  , data=cd[cd$year==2002,])
m.incumbency.04 <- felm(dv ~inc3 + ptynow + dpres + inc3*log(providers) + 
                          log(med.income) + log(pop) |state| 0 | state  , data=cd[cd$year==2004,])
m.incumbency.06 <- felm(dv ~inc3 + ptynow + dpres + inc3*log(providers) + 
                          log(med.income) + log(pop)|state| 0 | state  , data=cd[cd$year==2006,])
m.incumbency.08 <- felm(dv ~inc3 + ptynow + dpres + inc3*log(providers) + 
                          log(med.income) + log(pop)|state| 0 | state  , data=cd[cd$year==2008,])

#Creation of Table 10
stargazer(m.incumbency.02,m.incumbency.04,m.incumbency.06,m.incumbency.08, omit=c("cdfip", "year"),
          star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          covariate.labels = c("Incumbency$_{k}$", "Party Now$_{k}$","% Democratic Vote for President$_{k}$",
                               "log(Broadband Providers$_{k}$)","log(Median Income$_{k}$)",
                               "log(Population$_{k}$)", 
                               "Incumbency$_{kt}$*log(Broadband Providers$_{k}$)"),
          dep.var.labels   = "% Democratic Vote for House$_{kt}$",
          model.numbers = F,
          column.labels = c("2002", "2004", "2006", "2008"),
          add.lines = list(c("State F.E.", "Yes","Yes","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster(State)-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Incumbency Advantage",
          label="incumbency.cs.1")



m.pres.02 <- felm(inc.vote ~ factor(party)  +  log(providers)*pres.vote +
                    log(pop)  + log(med.income)|state| 0 | state , data=cd[cd$year==2002,])

m.pres.04 <- felm(inc.vote ~ factor(party)  +  log(providers)*pres.vote +
                    log(pop)  + log(med.income) |state| 0 | state  , data=cd[cd$year==2004,])

m.pres.06 <- felm(inc.vote ~ factor(party)  +  log(providers)*pres.vote +
                    log(pop)  + log(med.income) |state| 0 | state  , data=cd[cd$year==2006,])

m.pres.08 <- felm(inc.vote ~ factor(party)  +  log(providers)*pres.vote +
                    log(pop)  + log(med.income)  |state| 0 | state  , data=cd[cd$year==2008,])



stargazer(m.pres.02,m.pres.04,m.pres.06,m.pres.08, type="text")
#Creation of Table 11

stargazer(m.pres.02,m.pres.04,m.pres.06,m.pres.08, omit=c("cdfip", "year"),
          star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          covariate.labels = c("Republican$_{k}$", 
                               "log(Broadband Providers$_{k}$)",
                               "% Same Party Presidential Vote$_{k}$",
                               "log(Population$_{k}$)",
                               "log(Median Income$_{k}$)",
                               "% Same Party Presidential Vote$_{k}$*log(Broadband Providers$_{k}$)"),
          dep.var.labels   = "% Incumbent Vote$_{kt}$",
          model.numbers = F,
          column.labels = c("2002", "2004", "2006", "2008"),
          add.lines = list(c("State F.E.", "Yes","Yes","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster(State)-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Same Party Presidential Vote on Incumbent Vote Share",
          label="president.cs.1")



m.party.02 <- felm(inc.vote ~ factor(party)  +  log(providers)*party.unity +
                     log(pop)  + log(med.income)  |state| 0 | state , data=cd[cd$year==2002,])

m.party.04 <- felm(inc.vote ~ factor(party)  +  log(providers)*party.unity +
                     log(pop)  + log(med.income)  |state| 0 | state  , data=cd[cd$year==2004,])

m.party.06 <- felm(inc.vote ~ factor(party)  +  log(providers)*party.unity +
                     log(pop)  + log(med.income)|state| 0 | state  , data=cd[cd$year==2006,])

m.party.08 <- felm(inc.vote ~ factor(party)  +  log(providers)*party.unity +
                     log(pop)  + log(med.income) |state| 0 | state  , data=cd[cd$year==2008,])


#Creation of Table 12

stargazer(m.party.02,m.party.04,m.party.06,m.party.08, omit=c("cdfip", "year"),
          star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          covariate.labels = c("Republican$_{k}$", 
                               "log(Broadband Providers$_{k}$)",
                               "Party Unity$_{k}$",
                               "log(Population$_{k}$)",
                               "log(Median Income$_{k}$)",
                               "Party Unity$_{k}$*log(Broadband Providers$_{k}$)"),
          dep.var.labels   = "% Incumbent Vote$_{kt}$",
          model.numbers = F,
          column.labels = c("2002", "2004", "2006", "2008"),
          add.lines = list(c("State F.E.", "Yes","Yes","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster(State)-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Same Party Presidential Vote on Incumbent Vote Share",
          label="president.cs.1")


##### Appendix 5.3 Placebo Test ######
rm(list=ls())
load(file="PlaceboAggregateNatInc.Rdata")
load(file="MainModels.Rdata")



#cd[cd$year==2008 & cd$dv2==100,]



p.incumbency1 <- plm(dv ~ factor(year) + ptynow + dpres +inc3   ,
                     data=cd, model="within", index=c("cdfip","year"))
p.incumbency1.vcov <- vcovHC(p.incumbency1, type="HC0", cluster="group", method = "arellano")
p.incumbency1.se <- sqrt(diag(p.incumbency1.vcov))
p.incumbency2 <- plm(dv ~factor(year) + ptynow + dpres + inc3*log(providers)
                     + log(pop) + log(med.income),
                     data=cd, model="within", index=c("cdfip","year"))
p.incumbency2.vcov <- vcovHC(p.incumbency2, type="HC0", cluster="group", method = "arellano")
p.incumbency2.se <- sqrt(diag(p.incumbency2.vcov))


#Creation of Table 13

stargazer(p.incumbency1,p.incumbency2, omit=c("cdfip", "year"),
          se=list(p.incumbency1.se, p.incumbency2.se),star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          covariate.labels = c("Party Now$_{kt}$","% Democratic Vote for President$_{kt}$","Incumbency$_{kt}$",
                               "log(Broadband Providers$_{kt}$)","log(Median Income$_{kt}$)",
                               "log(Population$_{kt}$)", 
                               "Incumbency$_{kt}$*log(Broadband Providers$_{kt}$)"),
          dep.var.labels   = "% Democratic Vote for House$_{kt}$",
          model.numbers = F,
          add.lines = list(c("District F.E.", "Yes","Yes"),
                           c("Election Year F.E.","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Incumbency Advantage",
          label="placebo.incumbency")

#Marginal Effect Plot
set.seed(2105)
providers.sample  <- sample(cd$providers,500)
Rmarginsln <- function(model,xname,xzname,vcov.matrix,levels, ci=.95, latex=F){
  # Generate Marginal Effects
  m_effect <- coef(model)[xname]+coef(model)[xzname]*levels
  #Variances and Covariances
  v_x <- vcov.matrix[xname,xname]
  v_xz <- vcov.matrix[xzname,xzname]
  cov <- vcov.matrix[xname,xzname]
  #Standard Errors
  se <- sqrt(v_x + (levels^2)*v_xz+2*levels*cov)
  #T value for 95%CI
  if (ci==.95){
    t <- qt(0.025,model$df)
  }
  #T value for 90%CI
  if (ci==.9){
    t <- qt(0.05,model$df)
  }
  #Confidence Bounds
  m_upper <- m_effect + se*t
  m_lower <- m_effect - se*t
  # Remove Flotsom and Jetson
  #Printing Table
  table <- cbind(levels,m_effect,se,m_upper,m_lower)
  if (latex==T){
    library(xtable)
    print(xtable(table),include.rownames=FALSE)
  }
  if (latex==F){
    return(table)
  }
}


lnproviders <- log(seq(.1,20,.01))
margins.placebo <- Rmarginsln(p.incumbency2,"inc3", "inc3:log(providers)", p.incumbency2.vcov, levels = lnproviders)
margins.main <- Rmarginsln(m.incumbency2,"inc3", "inc3:log(providers)", m.incumbency2.vcov, levels = lnproviders)

#Creation of figure 2

pdf("./AppendixFigures/PlaceboIncumbencyAdvantageResults.pdf", height=6, width=8)
plot(exp(margins.placebo[,1]), margins.placebo[,2], type="n", ylim=c(3,12), xlim=c(2,16),
     xlab="Broadband Providers", ylab="Incumbency Advantage", lwd=2)
points(exp(margins.main[,1]), margins.main[,2], type="l", lwd=2, col="gray50")
lines(exp(margins.main[,1]), margins.main[,4], lty=2, col="gray50")
lines(exp(margins.main[,1]), margins.main[,5], lty=2, col="gray50")
points(exp(margins.placebo[,1]), margins.placebo[,2], type="l", lwd=2)
lines(exp(margins.placebo[,1]), margins.placebo[,4], lty=2)
lines(exp(margins.placebo[,1]), margins.placebo[,5], lty=2)
rug(providers.sample)
legend("topright", c("Main Effect","Placebo Effect"), lty=c(1,1), lwd=c(2,2), col=c("gray50","black"))
dev.off()


p.incwin.pres1 <- plm(inc.vote ~ factor(year)*factor(party) +   pres.vote
                      + log(pop) + log(med.income)
                      , data=cd, model="within", index=c("icpsr.id","year"))
p.incwin.pres1.vcov <-  vcovHC(p.incwin.pres1, type="HC0", cluster="group", method = "arellano")
p.incwin.pres1.se <- sqrt(diag(p.incwin.pres1.vcov))
summary(p.incwin.pres1)


p.incwin.pres2 <- plm(inc.vote ~ factor(year)*factor(party) +  log(providers)*pres.vote
                      + log(pop) + log(med.income)  , data=cd, model="within", index=c("icpsr.id","year"))
p.incwin.pres2.vcov <-  vcovHC(p.incwin.pres2, type="HC0", cluster="group", method = "arellano")
p.incwin.pres2.se <- sqrt(diag(p.incwin.pres2.vcov))
summary(p.incwin.pres2)


p.incwin.party1 <- plm(inc.vote ~ factor(year)*factor(party) +  party.unity 
                       + log(pop) + log(med.income), data=cd, model="within", index=c("icpsr.id","year"))
p.incwin.party1.vcov <-  vcovHC(p.incwin.party1, type="HC0", cluster="group", method = "arellano")
p.incwin.party1.se <- sqrt(diag(p.incwin.party1.vcov))


p.incwin.party2 <- plm(inc.vote ~factor(year)*factor(party) + log(providers)*party.unity
                       + log(pop) + log(med.income)   , data=cd, model="within", index=c("icpsr.id","year"))
p.incwin.party2.vcov <-  vcovHC(p.incwin.party2, type="HC0", cluster="group", method = "arellano")
p.incwin.party2.se <- sqrt(diag(p.incwin.party2.vcov))



#Creation of Table 14

stargazer(p.incwin.pres1, p.incwin.pres2,p.incwin.party1, p.incwin.party2, omit=c("icpsr.id","name","year","ptynow"),
          se=list(p.incwin.pres1.se, p.incwin.pres2.se,p.incwin.party1.se, p.incwin.party2.se),star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          order=c("pres.vote", "party.unity","log(providers)", "log(providers):pres.vote",
                  "log(providers):party.unity"),
          covariate.labels = c("Same Party Presidential Vote$_{jt}$",
                               "Same Party Presidential Vote$_{jt}$*log(Broadband Providers$_{jt}$)",
                               "Party Unity$_{jt}$",
                               "Party Unity$_{jt}$*log(Broadband Providers$_{jt}$)",
                               "log(Broadband Providers$_{jt}$)", 
                               "log(Population$_{jt}$)",
                               "log(Median Income$_{jt}$)" ),
          dep.var.labels   = "% Incumbent Vote$_{jt}$",
          model.numbers = F,
          add.lines = list(c("Member F.E.", "Yes","Yes","Yes","Yes"),
                           c("Election Year*Party F.E.","Yes","Yes","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Incumbency Advantage",
          label="incvote")






lnproviders <- log(seq(.1,20,.01))
margins.main <- Rmarginsln(m.incwin.pres2,"pres.vote", "log(providers):pres.vote", m.incwin.pres2.vcov, levels = lnproviders)
margins.placebo <- Rmarginsln(p.incwin.pres2,"pres.vote", "log(providers):pres.vote", p.incwin.pres2.vcov, levels = lnproviders)

#Creation of Figure 3a

pdf("./AppendixFigures/PlaceboPresVote.pdf", height=6, width=8)
plot(exp(margins.placebo[,1]), margins.placebo[,2], type="n", ylim=c(-.5,1), xlim=c(2,16),
     xlab="Broadband Providers", ylab="Effect of District Presidential Vote on Incumbent Margin", lwd=2)
points(exp(margins.main[,1]), margins.main[,2], type="l",col="gray50",lwd=2)
lines(exp(margins.main[,1]), margins.main[,4], lty=2,col="gray50")
lines(exp(margins.main[,1]), margins.main[,5], lty=2,col="gray50")
points(exp(margins.placebo[,1]), margins.placebo[,2], type="l",lwd=2)
lines(exp(margins.placebo[,1]), margins.placebo[,4], lty=2)
lines(exp(margins.placebo[,1]), margins.placebo[,5], lty=2)
rug(providers.sample)
legend("topright", c("Main Effect","Placebo Effect"), lty=c(1,1), lwd=c(2,2), col=c("gray50","black"))
#abline(h=0, lty=2)
dev.off()




lnproviders <- log(seq(.1,20,.01))
margins.main <- Rmarginsln(m.incwin.party2,"party.unity", "log(providers):party.unity", m.incwin.party2.vcov, levels = lnproviders)
margins.placebo <- Rmarginsln(p.incwin.party2,"party.unity", "log(providers):party.unity", p.incwin.party2.vcov, levels = lnproviders)

#Creation of Figure 3b

pdf("./AppendixFigures/PlaceboPartyUnity.pdf", height=6, width=8)
plot(exp(margins.placebo[,1]), margins.placebo[,2], type="n", ylim=c(-.6,.6), xlim=c(2,16),
     xlab="Broadband Providers", ylab="Effect of Party Unity on Incumbent Margin", lwd=2)
points(exp(margins.main[,1]), margins.main[,2], type="l",col="gray50",lwd=2)
lines(exp(margins.main[,1]), margins.main[,4], lty=2,col="gray50")
lines(exp(margins.main[,1]), margins.main[,5], lty=2,col="gray50")
points(exp(margins.placebo[,1]), margins.placebo[,2], type="l",lwd=2)
lines(exp(margins.placebo[,1]), margins.placebo[,4], lty=2)
lines(exp(margins.placebo[,1]), margins.placebo[,5], lty=2)
rug(providers.sample)
legend("topright", c("Main Effect","Placebo Effect"), lty=c(1,1), lwd=c(2,2), col=c("gray50","black"))
#abline(h=0, lty=2)
dev.off()

##### Appendix 6 Full Demographics Models####
rm(list=ls())
load(file="AggregateNatInc.Rdata")

m.incumbency1 <- plm(dv ~ factor(year) + ptynow + dpres + inc3 ,
                     data=cd, model="within", index=c("cdfip","year"))
m.incumbency1.vcov <- vcovHC(m.incumbency1, type="HC0", cluster="group", method = "arellano")
m.incumbency1.se <- sqrt(diag(m.incumbency1.vcov))
m.incumbency2 <- plm(dv ~factor(year) + ptynow + dpres + inc3*log(providers) + 
                       log(med.income) + log(pop) + med.age + perc.poverty + perc.white + perc.black + perc.bachelor ,
                     data=cd, model="within", index=c("cdfip","year"))
m.incumbency2.vcov <- vcovHC(m.incumbency2, type="HC0", cluster="group", method = "arellano")
m.incumbency2.se <- sqrt(diag(m.incumbency2.vcov))
summary(m.incumbency2)


#Creation of Table 15

stargazer(m.incumbency1,m.incumbency2, omit=c("cdfip", "year"),
          se=list(m.incumbency1.se, m.incumbency2.se),star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          covariate.labels = c("Party Now$_{kt}$","% Democratic Vote for President$_{kt}$","Incumbency$_{kt}$",
                               "log(Broadband Providers$_{kt}$)","log(Median Income$_{kt}$)",
                               "log(Population$_{kt}$)", "Median Age$_{kt}$",
                               "% Living in Poverty$_{kt}$","% White$_{kt}$",
                               "% Black$_{kt}$","% Bachelor Degrees$_{kt}$",
                               "Incumbency$_{kt}$*log(Broadband Providers$_{kt}$)"),
          dep.var.labels   = "% Democratic Vote for House$_{kt}$",
          model.numbers = F,
          add.lines = list(c("District F.E.", "Yes","Yes"),
                           c("Election Year F.E.","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Incumbency Advantage",
          label="incumbency.demo")







m.incwin.pres1 <- plm(inc.vote ~ factor(year)*factor(party) +   pres.vote +
                        + log(pop) + med.age + perc.poverty + perc.white + perc.black + perc.bachelor +  log(med.income)  , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.pres1.vcov <-  vcovHC(m.incwin.pres1, type="HC0", cluster="group", method = "arellano")
m.incwin.pres1.se <- sqrt(diag(m.incwin.pres1.vcov))
summary(m.incwin.pres1)


m.incwin.pres2 <- plm(inc.vote ~ factor(year)*factor(party) +  log(providers)*pres.vote +
                        + log(pop) + med.age + perc.poverty + perc.white + perc.black + perc.bachelor +  log(med.income)  , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.pres2.vcov <-  vcovHC(m.incwin.pres2, type="HC0", cluster="group", method = "arellano")
m.incwin.pres2.se <- sqrt(diag(m.incwin.pres2.vcov))
summary(m.incwin.pres2)


m.incwin.party1 <- plm(inc.vote ~ factor(year)*factor(party) +   party.unity + 
                         + log(pop) + med.age + perc.poverty + perc.white + perc.black + perc.bachelor +  log(med.income)  , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.party1.vcov <-  vcovHC(m.incwin.party1, type="HC0", cluster="group", method = "arellano")
m.incwin.party1.se <- sqrt(diag(m.incwin.party1.vcov))



m.incwin.party2 <- plm(inc.vote ~factor(year)*factor(party) + log(providers)*party.unity  + 
                         + log(pop) + med.age + perc.poverty + perc.white + perc.black + perc.bachelor +  log(med.income) , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.party2.vcov <-  vcovHC(m.incwin.party2, type="HC0", cluster="group", method = "arellano")
m.incwin.party2.se <- sqrt(diag(m.incwin.party2.vcov))





#Creation of Table 16

stargazer(m.incwin.pres1, m.incwin.pres2,m.incwin.party1, m.incwin.party2, omit=c("icpsr.id","name","year","ptynow"),
          se=list(m.incwin.pres1.se, m.incwin.pres2.se,m.incwin.party1.se, m.incwin.party2.se),star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          order=c("pres.vote", "party.unity","log(providers)", "log(providers):pres.vote",
                  "log(providers):party.unity"),
          covariate.labels = c("Same Party Presidential Vote$_{jt}$",
                               "Same Party Presidential Vote$_{jt}$*log(Broadband Providers$_{jt}$)",
                               "Party Unity$_{jt}$",
                               "Party Unity$_{jt}$*log(Broadband Providers$_{jt}$)",
                               "log(Broadband Providers$_{jt}$)", 
                               "log(Population$_{jt}$)",
                               "log(Median Income$_{jt}$)","Median Age$_{jt}$",
                               "% Living in Poverty$_{jt}$","% White$_{jt}$",
                               "% Black$_{jt}$","% Bachelor Degrees$_{jt}$"),
          dep.var.labels   = "% Incumbent Vote$_{jt}$",
          model.numbers = F,
          add.lines = list(c("Member F.E.", "Yes","Yes","Yes","Yes"),
                           c("Election Year*Party F.E.","Yes","Yes","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Incumbency Advantage",
          label="incvote.demo")


##### Appendix 7 No Demographics Models####
rm(list=ls())
load(file="AggregateNatInc.Rdata")

m.incumbency1 <- plm(dv ~ factor(year) + ptynow + dpres + inc3 ,
                     data=cd, model="within", index=c("cdfip","year"))
m.incumbency1.vcov <- vcovHC(m.incumbency1, type="HC0", cluster="group", method = "arellano")
m.incumbency1.se <- sqrt(diag(m.incumbency1.vcov))
m.incumbency2 <- plm(dv ~factor(year) + ptynow + dpres + inc3*log(providers),
                     data=cd, model="within", index=c("cdfip","year"))
m.incumbency2.vcov <- vcovHC(m.incumbency2, type="HC0", cluster="group", method = "arellano")
m.incumbency2.se <- sqrt(diag(m.incumbency2.vcov))
summary(m.incumbency2)

stargazer(m.incumbency1,m.incumbency2, omit=c("cdfip", "year"),
          se=list(m.incumbency1.se, m.incumbency2.se), type="text")

#Creation of Table 17

stargazer(m.incumbency1,m.incumbency2, omit=c("cdfip", "year"),
          se=list(m.incumbency1.se, m.incumbency2.se),star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          covariate.labels = c("Party Now$_{kt}$","% Democratic Vote for President$_{kt}$","Incumbency$_{kt}$",
                               "log(Broadband Providers$_{kt}$)",
                               "Incumbency$_{kt}$*log(Broadband Providers$_{kt}$)"),
          dep.var.labels   = "% Democratic Vote for House$_{kt}$",
          model.numbers = F,
          add.lines = list(c("District F.E.", "Yes","Yes"),
                           c("Election Year F.E.","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Incumbency Advantage",
          label="incumbency.demo")







m.incwin.pres1 <- plm(inc.vote ~ factor(year)*factor(party) +   pres.vote 
                      , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.pres1.vcov <-  vcovHC(m.incwin.pres1, type="HC0", cluster="group", method = "arellano")
m.incwin.pres1.se <- sqrt(diag(m.incwin.pres1.vcov))
summary(m.incwin.pres1)


m.incwin.pres2 <- plm(inc.vote ~ factor(year)*factor(party) +  log(providers)*pres.vote 
                      , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.pres2.vcov <-  vcovHC(m.incwin.pres2, type="HC0", cluster="group", method = "arellano")
m.incwin.pres2.se <- sqrt(diag(m.incwin.pres2.vcov))
summary(m.incwin.pres2)


m.incwin.party1 <- plm(inc.vote ~ factor(year)*factor(party) +   party.unity  
                       , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.party1.vcov <-  vcovHC(m.incwin.party1, type="HC0", cluster="group", method = "arellano")
m.incwin.party1.se <- sqrt(diag(m.incwin.party1.vcov))



m.incwin.party2 <- plm(inc.vote ~factor(year)*factor(party) + log(providers)*party.unity 
                       , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.party2.vcov <-  vcovHC(m.incwin.party2, type="HC0", cluster="group", method = "arellano")
m.incwin.party2.se <- sqrt(diag(m.incwin.party2.vcov))




stargazer(m.incwin.pres1, m.incwin.pres2,m.incwin.party1, m.incwin.party2, omit=c("icpsr.id","name","year","ptynow"),
          se=list(m.incwin.pres1.se, m.incwin.pres2.se,m.incwin.party1.se, m.incwin.party2.se), type="text")

#Creation of Table 18

stargazer(m.incwin.pres1, m.incwin.pres2,m.incwin.party1, m.incwin.party2, omit=c("icpsr.id","name","year","ptynow"),
          se=list(m.incwin.pres1.se, m.incwin.pres2.se,m.incwin.party1.se, m.incwin.party2.se),star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          order=c("pres.vote", "party.unity","log(providers)", "log(providers):pres.vote",
                  "log(providers):party.unity"),
          covariate.labels = c("Same Party Presidential Vote$_{jt}$",
                               "Same Party Presidential Vote$_{jt}$*log(Broadband Providers$_{jt}$)",
                               "Party Unity$_{jt}$",
                               "Party Unity$_{jt}$*log(Broadband Providers$_{jt}$)" ),
          dep.var.labels   = "% Incumbent Vote$_{jt}$",
          model.numbers = F,
          add.lines = list(c("Member F.E.", "Yes","Yes","Yes","Yes"),
                           c("Election Year*Party F.E.","Yes","Yes","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Incumbency Advantage",
          label="incvote.demo")


